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Abstract. Aftershock sequences are of particular interest in seismic re- 
search since they may condition seismic activity in a given region over long 
time spans. While they are typically identified with periods of enhanced seis- 
mic activity after a large earthquake as characterized by the Omori law, our 
knowledge of the spatiotemporal correlations between events in an aftershock 
sequence is limited. Here, we study the spatiotemporal correlations of two 
aftershock sequences form California (Parkfield and Hector Mine) using the 
recently introduced concept of "recurrent" events. We find that both sequences 
have very similar properties and that most of them are captured by the space- 
time epidemic-type aftershock sequence (ETAS) model if one takes into ac- 
count catalog incompleteness. However, the stochastic model does not cap- 
ture the spatiotemporal correlations leading to the observed structure of seis- 
micity on small spatial scales. 
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1. Introduction 

One of the grand challenges for seismology is to establish the relationship between 
stress and strain in the lithosphere [Forsyth et al, 2009]. While motions of tectonic plates 
and surface deformations can be measured precisely with satellite imaging and networks 
of Global Positioning System receivers, strainmeters, seismometers and tiltmeters, the 
causative stress can only be inferred. However, the temporally and spatially dependent 
rheology which describes the linkage between the forces (stresses) and the resulting defor- 
mations (strains) is generally not well understood [Kanamori and Brodsky, 2004]. This 
makes it currently impossible to conclusively answer how some earthquakes trigger other 
earthquakes thousands of kilometers away, for example, and to predict earthquakes reli- 
ably. 

An alternative approach to gain insight into the underlying dynamics of earthquakes 
is to study the spatiotemporal patterns of seismicity where each earthquakes is treated 
as a point event in space and time with a given magnitude. Such an approach may 
shed light on the fundamental physics since these patterns are emergent processes of 
the underlying many-body nonlinear system. Indeed, it has been proved successful in 
many cases and led to the discovery of new key features of seismicity [Rundle et al, 
2003; Corral, 2004; Davidsen and Goltz, 2004; Shcherbakov et al, 2004; Davidsen and 
Paczuski, 2005; Baiesi and Paczuski, 2005; Felzer and Brodsky, 2006; Hainzl et al, 2006; 
Corral, 2006; Marsan and Lengline, 2008; Zaliapin et al, 2008]. For example, in Davidsen 
et al. [2006, 2008] the spatiotemporal clustering of earthquakes in Southern California was 
found to show non-trivial features which led to a new and independent estimate of the 
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rupture length and its scaling with magnitude. In particular, the results provide further 

evidence for a shadowing effect associated with smaller earthquakes [Rubin, 2002; Fischer 
and Hordlek, 2005; Hainzl and Marsan, 2008]. The key to these findings was a unique 
approach to quantify non-trivial spatiotemporal clustering based on the view that any 
suitable definition of clustering should be purely contextual and depend only on the actual 
history of events without any further assumptions [Davidsen et ai, 2008]. This approach 
utilizes the notion of space-time records to define "recurrences" and maps seismicity onto 
a graph or network, thus, allowing the characterization of spatiotemporal clustering by 
means of tools from complex network theory [Albert and Barabasi, 2002; Newman, 2003]. 

To elucidate the origin of the observed non-trivial clustering further, we study here 
the spatiotemporal correlations of aftershock sequences which follow large shallow earth- 
quakes. Aftershocks are the most obvious example of earthquakes that are triggered in 
part by preceding events as follows from the observed increased seismic activity captured 
by the Omori law [Utsu et ai, 1995]. Thus, their specific clustering in space and time 
should provide information on the underlying dynamics and triggering mechanisms — see, 
for example, [Felzer and Brodsky, 2006; Gomberg and Felzer, 2008]. Moreover, aftershocks 
are important from a conceptual point of view since the current main paradigm in statisti- 
cal seismology classifies earthquakes as triggered events like aftershocks and "background" 
events which are hypothesized to be induced by other means. It is important to realize, 
however, that the notion of background events — and aftershocks, for that matter — is 
not well-defined as no clear physical difference between such events has been established. 
In particular, background events could be artifacts to a large extent since small earth- 
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quakes can trigger larger events [Helmstetter et ai, 2006; Marsan and Lengline, 2008] and 

many small earthquakes are typically not detected [Sornette and Werner, 2005a, b]. 

For the aftershock sequences associated with the Parkfield earthquake (September 28, 

2004) and the Hector Mine earthquake (October 16, 1999), we find that both sequences 

show very similar spatiotemporal correlations as quantified by the method of Davidsen 

et al. [2006, 2008]. While most features are also similar to what has been observed 

previously for a 20- year catalog from Southern California [Davidsen et ai, 2008], there 

is no indication of a shadowing effect. Moreover, we find that most of the observed 

features can be captured by a space-time version of the epidemic-type aftershock sequence 

(ETAS) model [Helmstetter and Sornette, 2002a; Ogata and Zhuang, 2006] if catalog 

incompleteness [Lennartz et al, 2008a] is taken into account. This suggests that the 

spatiotemporal correlations are to a large extent a consequence of a few established laws 

of seismicity and missing data. This was further confirmed by a comparison with a 

simple non-homogeneous Poisson (NHP) model following the Omori law, but without any 

spatiotemporal correlations between events. Yet, we find that this stochastic model has 

several shortcomings with respect to the spatial distribution of seismicity. 

2. Aftershock sequences and the ETAS model 

Although there is considerable statistical variability associated with aftershocks, their 
behavior appears to satisfy several scaling laws to a reasonably good approximation. 
Among them are the Gutenberg-Richter scaling relation for the frequency-magnitude 
distribution [Gutenberg and Richter, 1949] which is certainly satisfied on long time 
scales [Shcherbakov et al, 2006], Bath's law for the difference between the magnitudes 
of a mainshock and its largest aftershock [Bath, 1965] as well as the modified Omori law 
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for the temporal decay of aftershock rates [Utsu et al, 1995; Shcherbakov et al, 2004; 

Hainzl and Marsan, 2008; Narteau et al, 2009]. More recently, another scaling law char- 
acterizing the epicenter distribution of aftershocks has been found [Felzer and Brodsky, 
2006]. 

In an attempt to establish a statistical null model of aftershocks which would incorporate 
some of these scaling laws, the epidemic-type aftershock sequence (ETAS) model was 
introduced [Kagan and Knopoff , 1987; Ogata, 1988; Helmstetter and Sornette, 2002b]. 
This model describes a stochastic branching process in which any earthquake may trigger 
other earthquakes, which in turn may trigger more, and so on. In particular, the location 
and the time of occurrence of each "daughter" event is strongly correlated to its "mother" 
event: The occurrence rate of daughter events at time t and position r triggered by a 
mother event of magnitude m,, at time ti and position r iy is defined as 

(PmXt-t^r-fi) = p(mi)ty(t - ti)<f> mi (r - fj), (1) 

where p{mi) is the average number of aftershocks directly triggered by an event of mag- 
nitude rrii, ty(t — t^ is the normalized temporal distribution of directly triggered events 
and $ mi (r — r*j) is the normalized spatial distribution of directly triggered events in two 
dimensions. p{mi) is assumed to follow the productivity law 

p{mi) = Kl0 a{m >~ mo \ (2) 

where m 8 > m and K, a are constants. Note that the productivity law is zero below m 
implying that events smaller than m do not trigger other earthquakes. This condition is 
necessary to ensure a finite total occurrence rate — see [Sornette and Werner, 2005a] for 
a discussion of its physical relevance. ty(t) is assumed to be determined by the modified 
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Omori law which can be written as 



0c e 

(t + c y 



*(*) = u , ( 3 ) 



where c, # > and is the Heaviside step function. Note that the exponent 1 + 9, 
which describes the time distribution of the direct aftershocks, is typically larger than 
the observed Omori exponent, which characterizes the whole cascade of directly and in- 
directly triggered aftershocks [Helmstetter and Sornette, 2002b; Marsan and Lengline, 
2008]. Finally, $ m (r) is assumed to follow the recently established epicenter distribution 
of aftershocks [Helmstetter and Sornette, 2002b; Felzer and Brodsky, 2006] which can be 
expressed as 

$m(r) = d(m)(\r\/d(m) + l) 1+ ^ (4) 
where \i = 0.35 and d{m) = <i 10°' 45m with d = 15m is the rupture length of the triggering 

event of magnitude m [Wells and Coppersmith, 1994; Kagan, 2002; Davidsen et ai, 2008]. 

Note that this expression is already problematic since it assumes a rotationally symmetric 

distribution of aftershocks, thus, neglecting the typically anisotropic fault structure on 

which aftershocks occur and also neglecting the frequent observation that the mother 

event is located on the margin of the area of triggered events. Since the epicenters of the 

aftershock sequences considered here are approximately located on a single fault, we focus 

on the one-dimensional case and discuss the two-dimensional case as necessary. 

The magnitude of each event is independently sampled from the normalized 

Gutenberg-Richter distribution, 

P{m) = Mn(10)10- 6(m - mo) , (5) 
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where the constant b is typically close to 1 and m is again the lower threshold, below 

which no aftershocks are initiated. 

An important quantity in the ETAS model is the average number of n of daughter 

events per earthquake, averaged over all magnitudes, which is given by 



where the last equation assumes a < b which has been suggested as the relevant 
regime [Marsan and Lengline, 2008]. 

2.1. Data sets and parameter values 

Here, we study the aftershock sequences associated with the Parkfield earthquake 
(M = 6.0, September 28, 2004) and the Hector Mine earthquake (M = 7.1, October 16, 
1999) as identified in Shcherbakov et al [2004, 2006]. In both cases, the high-quality seis- 
mic network in the vicinity of these mainshocks provided a particularly well-documented 
sequence of aftershocks. This is especially true for the Parkfield sequence [Bakun et al, 
2005]. For Parkfield, the number of identified aftershock over a time period of T = 365 
days is 2056 above magnitude m = 1.15. For Hector mine, we have 5380 aftershocks above 
m = 2.0. Both aftershock sequences differ significantly in their spatial extent and also 
slightly in the exponents of the Omori law and the Gutenberg-Richter law. 

In principle, it seems straightforward to estimate the parameters of the ETAS model 
for a given aftershock sequence. Yet, an important aspect of any aftershock sequence is 
that at early times after a mainshock not all (small) aftershocks are detected due to the 
large amount of seismic noise [Kagan, 2004; Kagan and Houston, 2005; Peng et al, 2006]. 
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This has led to the proposition of a time-dependent magnitude threshold of completeness 

which involves further parameters [Helmstetter et al, 2006]. To take this into account, we 
proceed as follows: First we generate artificial aftershock sequence from the ETAS model 
for a given main event with magnitude M as outlined in the Appendix Al. Then, we 
impose two constraints on the generated catalog to take into account the finite area size 
covered by the empirical catalogs and to mimic missing data: i) events that lie outside of 
the spatial area considered for the empirical catalogs were discarded; ii) events below the 
time- varying magnitude threshold of completeness m c (M, t) were discarded with a certain 
probability according to the procedure described in Lennartz et al. [2008b]. To be more 
precise, we define mc(M, t) = max [mi(M, t), m 2 ] where mi (M, t) = M— 4.5 — 0.75 log 10 (t) 
reflects the incompleteness at time t after the main event of magnitude M, and the time- 
independent sensitivity of the seismic network, different for each empirical catalog, is 
captured by the constant m 2 . Then, the probability of observing an event of magnitude 
smaller than the threshold m c (M,t) is given by Pj{m) 

= iQ-iMi-™)^ w here j = 1 if 

m c (M,t) = m 1 (M,t) and j = 2 otherwise. 

Using the parameters estimated in Lennartz et al. [2008a] for the Parkfield and Hector 
mine aftershock sequences, only K and a of the ETAS model as well as the fraction of 
background events have to be approximated. Our selection criteria were both the total 
number of events in the catalogs, and the functional form of the total event rate (see 
Eq. (Al)). The background rate was set at 1 event /day, which gives a total event rate in 
very good agreement with the empirical data (not shown). The other two parameters were 
chosen, such that the average catalog size was similar to the empirical catalogs (after some 
events were discarded by the procedure outlined above), i.e., 2056 for Parkfield and 5380 
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for Hector Mine, with an accepted deviation of at most ±200 events per catalog. In total, 

100 artificial Parkfield and Hector Mine catalogs were generated. The parameter values 

used for simulating the two different aftershock sequences are summarized in table 1. 

3. Network of recurrences 

We analyze the two aftershock sequences and the catalogs generated by the ETAS 
model according to a recently proposed method which allows one to characterize the 
spatiotemporal clustering of seismicity [Davidsen et ai, 2006, 2008]. It is based on the 
notion of a recurrence in the context of spatiotemporal point processes: An event is defined 
to be a recurrence of a given previous event if it occurs closer in space to that event than 
all intervening events. Recurrences are therefore record breaking events with respect to 
distance. This relationship is naturally represented by a directed network of recurrences: 
Each event a k defined by its location r k and time of occurrence t k , with k — 1, ...,N, is 
a vertex in the network and a directed edge from a k to a m exists for k < m if 
recurrence of a k , i.e., if the distance \ f m — r k \ is smaller than the distance \r k > — r k \ for all 
events a k > with k < k! < m. This definition assumes that the events are ordered according 
to their occurrence in time. Each recurrence, i.e. each edge on the network, is therefore 
characterized by the time interval t — t m — t k between the two connected events k and 
m and analogously by the spatial distance / between the two. Note that the mapping of 
the point process dynamics to the recurrence network is entirely data-driven and does not 
impose any arbitrary space and time scales other than those associated with the given 
event catalog. 
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To investigate the dependence of the network properties on the magnitude thresholds of 

events considered, we analyze networks obtained for different magnitude thresholds. This 

is crucial to identify robust features as well as potential scaling properties. 

4. Results 

4.1. Topological characteristics of the network 

We turn now to the analysis of the statistical properties of the network of recurrences 
which have been proved useful in the analysis of the overall seismic activity in Southern 
California [Davidsen et al, 2006, 2008]. The growth of the network can be measured 
by the average degree (k) as a function of the number of events N in the catalog. The 
in(out)-degree of a node is the number of edges pointing towards (away from) it and the 
in- and out-degree averaged over the entire network are obviously the same. The number 
of nodes iV can be controlled by changing the magnitude threshold above which events 
are considered for constructing the network. As can be seen in Fig. 1, the ETAS datasets 
as well as the empirical catalogs show a logarithmic increase of (k) with N. Both for 
Parkfield and Hector Mine, the increase is compatible with the ETAS data sets if one 
takes into account the statistical uncertainties. In all cases, the behaviour is similar to 
what is expected if events were randomly, independently and uniformly placed in space and 
time, which would imply that the growth should follow (k) ~ In N for large N [Davidsen 
et al, 2008]. Actually, this result even holds if the events are not uniformly distributed in 
time including the case of a simple non-homogeneous Poisson (NHP) model following the 
Omori law which by construction does not have any spatiotemporal correlations between 
events. The comparison with such a model is particularly helpful to identify statistical 
properties of the empirical aftershock series and of series generated by the ETAS model 
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that arise trivially and are, thus, not due to spatiotemporal correlations. The properties 

of the NHP model are analytically derived and discussed in the Appendix A2. 

The probability distributions of in- and out-degrees can be seen in Fig 2. The out-degree 
distributions for both empirical and ETAS catalogs are in excellent agreement. While the 
distributions are compatible with a Poisson distribution over the range where we have 
non-zero values for the respective empirical catalog, the decay in the distribution for the 
respective ETAS catalog for larger values of k out is slightly but significantly slower than a 
Poisson distribution (not shown). The reason that we actually observe larger values of k out 
in the ETAS catalogs is due to the much better statistics. The situation is different for the 
in-degree distributions. Fig 2 shows that there are deviations between the distributions 
for the empirical catalog and for the ETAS catalog for both Parkfield and Hector mine. 
The deviations are most pronounced for large k in and small m. In comparison to a Poisson 
distribution, the distributions are more narrow (not shown). 

For the NHP model, a Poisson distribution is expected both for the in- and out-degree 
distributions in the limit of large networks. The observed deviations imply that the de- 
gree distributions capture some of the underlying non-trivial spatiotemporal correlations. 
Moreover, the differences in the in-degree distribution between the empirical data and the 
ETAS model are a first indication that the spatiotemporal correlations generated by the 
ETAS model deviate from the empirical ones. 

4.2. Temporal recurrence statistics 

We turn now to the statistics of recurrence time intervals t, which are associated with 
the edges of the recurrence network. The comparison between the empirical catalogs and 
the ETAS model is shown in Fig 3. It can be seen that, with the parameters used for the 



DRAFT 



April 13, 2010, 1:29am 



DRAFT 



PEIXOTO ET AL.: SPATIOTEMPORAL CORRELATIONS OF AFTERSHOCK SEQUENCES X - 13 

ETAS model, and accounting properly for the missing data, the agreement between the 
datasets is excellent. Note that deviations for short time intervals are largely within the 
statistical uncertainties, but a lack of resolution and missing data in the empirical data 
not accounted for certainly play a role as well. This is further supported by the fact that 
there are basically no significant deviations for the Parkfield sequence which benefits from 
the higher quality of the seismic network in the region. 

As the insets of Fig 3 indicate, after a characteristic time, which is independent of the 
magnitude threshold m, the probability density functions (PDF) p(t) of the time intervals 
decay as t~ a9 and t~ 1A for Parkfield and Hector Mine, respectively. If one does not account 
for missing data in the ETAS models, the exponents of the power-law decay are unchanged 
but there are significant and systematic deviations at small time scales (not shown). In 
particular, the characteristic time varies with m. This is not only further evidence that 
it is crucial to take the effect of missing data into account but it also indicates that the 
characteristic time — which is of the order of a minute — might be an observational 
artifact. 

Further support comes from the NHP model which predicts p(t) as shown in the upper 
panel of Fig. 3 for Parkfield and different magnitude thresholds m (see Appendix A2 for 
details). Indeed, the apparent lack of variation in p(t) with respect to m for the NHP 
model is due to the dependence of the parameters in the Omori law given in Eq. (A2) on 
the magnitude threshold. In particular, the change in c — which is typically related to 
missing data — prevents any significant variation in the characteristic time which would 
be expected otherwise (see Eq. (A9)). However, the NHP model also shows a faster decay 
than what is actually observed for larger values of t. In fact, Eq. (A9) predicts that we 
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should not expect a pure power-law for Parkfield but rather a power-law with logarithmic 

corrections. We attribute this difference between the NHP model on one side and the 

ETAS model and the empirical data on the other side to the presence of spatiotemporal 

correlations. Clearly, this implies that at least some of the spatiotemporal correlations 

present in the empirical data are correctly captured by the ETAS model. 

4.3. Spatial recurrence statistics 

In addition to the time interval, there is also the spatial distance I associated with 
each edge of the network. The PDF of these distances p(l) for both the empirical and 
ETAS catalogs are shown in Figs. 4 and 5 for different magnitude thresholds m and 
different numbers of events N, respectively. Contrary to the temporal statistics, there are 
significant differences between the ETAS model and the empirical data. 

For the empirical data, p(l) increases for small distances up to a characteristic distance 
followed by a decay for larger arguments 1 . The decaying part itself consists of two 
regimes, a power-law decay for intermediate and large distances and a much faster decay 
for the largest arguments. The latter is simply a finite size effect since the maximal spatial 
extent of the aftershock sequences is about 40km for Parkfield and about 150km for Hector 
Mine. For both Parkfield and Hector mine, the power-law decay has the same exponent 
of about 1.05 as evident from the insets of Figs. 4 and 5. The inset of Fig. 4 also shows 
that in both cases the characteristic distance scales as 10 45m with magnitude threshold 
m, while it scales as N~ 0A5 with the number of considered events iV (inset of Fig. 5). 
Since iV oc 10~ m according to the Gutenberg- Richter law for b = 1, both scaling laws 
are trivially related in this case. As a comparison of the insets of Fig. 5 shows, the only 
obvious difference between the aftershock sequences of Parkfield and Hector mine are the 
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constant prefactors in the respective scaling laws. This can be attributed to the different 

geometry of the faults in the areas considered. It is important to realize that all these 

findings are not significantly different from those for the NHP model (see Eq. (A7, A8)) 

suggesting that they do not reflect any non-trivial spatiotemporal correlations. 

For the ETAS model (see Figs. 4 and 5), p(l) shows a similar behavior for large and 

intermediate distances but there are huge deviations at length scales smaller than the 

characteristic distance discussed above. While there is a characteristic distance for the 

ETAS model, it scales as A -0 - 8 and 10 a8m , respectively 2 . This is not only different from 

the empirical data but there is also no obvious relationship to the spatial propagation of 

activity in the ETAS model defined by Eq. (4). Moreover, p{l) does not increase for small 

distances but is clearly constant. Some of these differences can be reasonably attributed 

to the different dimensionality of the datasets: We have assumed for the ETAS model 

that the epicenters of the events are displaced in one dimension, in order to mimic the 

displacement along a single fault. However, this assumption breaks down exactly at small 

distances, since on those scales the higher-dimensional structure of the fault itself becomes 

important. To address this point, we have analyzed a version of the ETAS model for 

Parkfield where the events are distributed in a two-dimensional area. While this change 

in dimensionality does not affect the overall functional form of p(l), the characteristic 

distance now scales as N~ 0A5 and l0°- 45m , respectively. Despite this agreement with the 

empirical data, the prefactor in the scaling of the characteristic distance is more than an 

order of magnitude larger. It is also important to note that p(l) remains constant for 

small distances. This indicates that the spatial version of the ETAS model considered 

here is not able to reliably recover the spatial distribution of aftershocks on small scales. 
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This is not unexpected since the isotropic spatial distribution of triggered events used in 

the ETAS model does not account for the orientations of the different faults and rupture 

areas. Surprisingly, however, the functional form of p(l) and the scaling behavior of the 

characteristic distance observed for the empirical data can be well-reproduced by the 

NHP model (see Eq. (A7, A8)) despite an isotropic spatial distribution of events in two 

dimensions. This is an indication that not necessarily the isotropy of the normalized 

spatial distribution of directly triggered events given in Eq. (4) is the issue but rather 

those spatiotemporal correlations that affect short spatial scales. 

5. Discussion 

It is interesting to compare the properties of the network of recurrences for the after- 
shock sequences associated with the Parkfield earthquake (September 28, 2004) and the 
Hector Mine earthquake (October 16, 1999) with those found for a 20-year catalog from 
southern California [Davidsen et al, 2008]. For example, the out-degree distributions are 
clearly different: While they seem to follow a Poisson distribution for the isolated after- 
shock sequences, the distributions are significantly broader for the 20-year catalog. This 
indicates that spatiotemporal correlations exist between earthquakes that extend beyond 
what is typically considered an aftershock sequence. Another important finding is that 
the distributions of the time intervals between events and their recurrences are basically 
indistinguishable. Since the analysis of the ETAS model indicates that the absence of any 
dependence on the threshold magnitude for the aftershock sequences is a consequence of 
missing data, one can draw an analogous conclusion for the 20-year catalog. Finally, it is 
important to realize that the distributions of the spatial distance between an earthquake 
and its recurrences for the aftershock sequences do not show any indication of a non-trivial 
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scaling with threshold magnitude. This is in sharp contrast to the 20-year catalog, which 

provided further evidence for a shadowing effect associated with smaller earthquakes [Ru- 
bin, 2002; Fischer and Hordlek, 2005; Hainzl and Marsan, 2008]. However, the absence of 
a shadowing effect might be attributed to the fact that the main events of both aftershock 
sequences studied here were of large magnitude — for which stress shadows are seldom 
observed [Helmstetter et al, 2005; Main, 2006] - and, hence, dominated the statistical 
behavior. 

The comparison of the non-trivial statistical features of the network of recurrences for 
the aftershock sequences with those for the ETAS model showed that most of them are 
indeed correctly captured by the model if the effect of catalog incompleteness is taken into 
account. While this is a clear indication that some of the spatiotemporal correlations be- 
tween earthquakes in aftershock sequences can be well-understood within the framework 
of the stochastic model, the significant differences in the spatial statistics of recurrences 
imply that the ETAS model is inadequate to explain all spatiotemporal correlations associ- 
ated with aftershock sequences. Thus, our analysis confirms that the statistical properties 
of the network of recurrences can serve as a valuable benchmark test for models of seis- 
micity as found in an earlier work on a conceptual model of earthquake dynamics [Peixoto 
and Davidsen, 2008]. It remains to be seen whether other stochastic models of seismicity 
pass this test [Vere- Jones, 2005; Turcotte et al, 2007]. 

Appendix A 

Al. Numerical simulation of the ETAS model 

An artificial earthquake catalog can be created in an efficient manner with the ETAS 
model as follows. Starting from an initial event of magnitude M at the origin, and at 
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time t = 0, the number of daughter events is sampled from a Poisson distribution with 

average p(M). For each daughter event i, its magnitude, position, and occurrence time 

are independently sampled from P(mj), $(fj,mj) and ^(tj — t), respectively. The same 

procedure is then repeated for each daughter event, recursively, until there are no more 

events to be generated up to a predetermined maximum time. Background events, which 

are not daughters of any other event, can also be added, and their daughters can be 

obtained in the same manner. Thus, the total number of iterations is proportional to 

the number of events in the final catalog. The final rate of occurrence at a given time, 

considering all daughter events, is then given by 



where K = n(b — a)/b, \ b is the constant rate of background events, and the sum is taken 
over all previous events in the catalog. 

A2. Simple non-homogeneous Poisson model 

We follow the same approach as in Davidsen et al. [2008] which allows us to define a 
simple non-homogeneous Poisson (NHP) model by specifying the functions /!;(/) and Ht(t). 
These functions quantify the spatial and temporal distribution of events with respect to 
some reference event, respectively. Note that this fact trivially excludes the existence of 
any spatiotemporal correlations between events. In contrast to [Davidsen et al, 2008], we 
assume here that the rate of activity is not constant but given by the modified Omori law, 



where we restrict ourselves to the case p > 1 as observed for Parkfield [Shcherbakov et al, 
2006]. We note that p is not the same as 9 + 1 in Eq. (3) for the ETAS model, since the 
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latter describes the rate of events directly triggered by a single event, while the former 

describes the total event rate. Despite this non-homogeneous rate of activity, some of 

the statistical properties of the network of recurrences remain unchanged, including the 

topological structure, since they only depend on the spatial characteristics of the point 

process and the temporal ordering [Davidsen et ai, 2008]. Incorporating further the effect 

of finite space-time domains, we have 



fjn(l) = 2alB (L-l), 



K 



Q( T _t-t ). 



(A3) 
(A4) 



{c + t + to) 1 

Here, a is some constant such that /i/(7) describes a homogeneous distribution of events 
in a two-dimensional ball of radius L around the reference event. Hence, we assume 
translational invariance in space. T corresponds to the finite observation period after the 
mainshock occurring at time t = 0. to is the time of occurrence of the reference event. 
Using the formalism presented in Davidsen et al. [2008], one can compute the average 
PDFs of finding a recurrence at distance / from or at time t after a randomly chosen 
reference event 3 . This leads to 
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The functional behavior of these PDFs can be broken down into several parts. For pi(l), 

we find 



aK 



Pi(0 



' (p-i) \ cP 
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(c+T) 



I* < I < L 
I > L, 



(A7) 



where the characteristic distance /* scales as 

1 ( 1 
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2(p-l)\cP- 1 (c + T) p -\ 
Note that the total number of events N is simply given by A^ = aL 2 J Q T n(t)dt. A com- 
parison with the results derived in Ref. [Davidsen et al, 2008] for a homogeneous rate of 
activity shows that the qualitative features and the scaling exponents are unchanged. For 
Pt(t), we have 
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see Eq. (A10) 
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where r(-,-) denotes the incomplete Gamma function 4 . Note that depending on the 
exact parameter values the 1/t regime might be present or not 5 . For the aftershock 
sequence of Parkfield, one can also obtain p t (t) by numerically integrating Eq. (A6). The 
values of the constants in Eqs. (A3,A4) can be derived from the values estimated for 
the Omori law given in Shcherbakov et al. [2006]. Note that the values of c and K vary 
with the magnitude threshold m. The variation in c can arise due to missing data - 
see Shcherbakov et al. [2004, 2006], however, for a different interpretation. For m = 2.15, 
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we have c = 0.016 d and aL 2 K = 21 d p_1 . For m = 1.15, we have c = 0.08 d and 

aL 2 K = 200 d p_1 . These curves are plotted in Fig. 3. 

For T — > oo and L — > oo, Eqs. (A6,A9) reduce to 

f°° p — 1 

Mt) = <p- i) <f l J t {c + toHc + t + tor _ {c + tof{c + t + to) ^ (All) 

p,(f) ~{h^[m© + ^] <>>r, (A12) 

where t* oc c. 
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Notes 

1. The relative location errors of the epicenters are typically between 10's of meters and 100's of meters for both aftershock 
sequences. 

2. We note that in Figs 4 and 5 distances smaller than 10 km were not considered, in order to obtain a better comparison 
with the empirical datasets. However, in the case of the ETAS model, this region is statistically significant, and removing 
it from the distribution destroys any chance of a symmetric collapse by a uniform scaling of both x and y axes. In our 
case, it suffices to scale both axes independently, where the scaling of the rr-axis gives the characteristic length exponent, 
and the j/-axis scaling is done purely in an ad hoc manner so that the curves collapse. 

3. It is important to realize that the average has to take into account that the rate of activity decays with time 

4. The above approximation for t 2> c is based on splitting the integral into two parts, f °° . . . dto = /o C • ■ ■ + / c °f t 

The first term was then approximated assuming (c + to) <S t and taking only terms of first order for (c + to) into account. 

The second term was approximated assuming (c + to) ^> t, again taking only first order terms of t into account. 

5. There are more cases that can be considered for pt(t), especially one for t > (T + c)/2, which describes a sharp cutoff near 

t. This case, however, only occurs over one order of magnitude and is therefore not important for a scaling-comparison 
with empirical data. 
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Event M m Q 9 c (days) b a K ji m 2 72 (days x ) 

Parkfield 6.0 1.15 0.09 0.00395 0.89 0.8 0.2 0.7 1.2 1.5 1 

Hector Mine 7.1 2.0 0.21 0.024 1.01 0.789 0.28 1.8 1.4 1.5 1 
Table 1. Parameters used for the ETAS model. 
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Figure 1. Average degree (k) as a function of the number of events N used to build 
the network. The values of iV were obtained by increasing the magnitude threshold m of 
the event selection. The label above each plot indicate the corresponding m values for 
the empirical data sets only. The error bars correspond to the standard deviation for the 
empirical datasets, and to the standard deviation of the average for the ETAS datasets. 
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Figure 2. In- and out-degree distributions for different magnitude thresholds m, for 
Parkfield (above) and Hector Mine (below). Solid lines represent results obtained with 
the ETAS model, and symbols the empirical catalogs. 
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Figure 3. Probability density functions of the time intervals between an event and its 
recurrences for different threshold magnitudes m, for Parkfield (above) and Hector Mine 
(below). Solid lines represent results obtained for the ETAS model, and symbols for the 
empirical catalogs. The additional solid lines in the upper panel correspond to curves 
obtained with the NHP model. The insets show the rescaled functions as indicated by 
the axis labels. 
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Figure 4. Recurrence distance distribution for different threshold magnitudes m, for 
Parkfield (above) and Hector Mine (below). Solid lines represent the respective ETAS 
model, and symbols correspond to the empirical catalogs. The insets show the rescaled 
distributions, with j x — ly — 0.45 for both empirical data sets, 7 X = 0.8 for both ETAS 
data sets, 7 y = —0.1 for the ETAS model of the Parkfield sequence and 7 y = —0.15 for 
the ETAS model of the Hector Mine sequence 1 . 
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Figure 5. Recurrence distance distributions if only the first N events after the main- 
shock are considered, for Parkfield with threshold magnitude m = 1.15 (above) and Hector 
Mine with m = 2 (below). Solid lines represent the respective ETAS model, and symbols 
correspond to the empirical catalogs. The insets show the rescaled distributions, with 
7x — ly — 0.45 for both empirical data sets, 7 X = 0.8 for both ETAS data sets, 7 y = —0.1 
for the ETAS model of the Parkfield sequence and jy = —0.15 for the ETAS model of the 
Hector Mine sequence 1 . 



DRAFT 



April 13, 2010, 1:29am 



DRAFT 



